home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / DJLSR106.ARJ / ACG.CC < prev    next >
C/C++ Source or Header  |  1992-03-30  |  10KB  |  292 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /* 
  3. Copyright (C) 1989 Free Software Foundation
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17.  
  18. #ifdef __GNUG__
  19. #pragma implementation "ACG.h"
  20. #endif
  21. #include <ACG.h>
  22. #include <assert.h>
  23.  
  24. //
  25. //    This is an extension of the older implementation of Algorithm M
  26. //    which I previously supplied. The main difference between this
  27. //    version and the old code are:
  28. //
  29. //        + Andres searched high & low for good constants for
  30. //          the LCG.
  31. //
  32. //        + theres more bit chopping going on.
  33. //
  34. //    The following contains his comments.
  35. //
  36. //    agn@UNH.CS.CMU.EDU sez..
  37. //    
  38. //    The generator below is based on 2 well known
  39. //    methods: Linear Congruential (LCGs) and Additive
  40. //    Congruential generators (ACGs).
  41. //    
  42. //    The LCG produces the longest possible sequence
  43. //    of 32 bit random numbers, each being unique in
  44. //    that sequence (it has only 32 bits of state).
  45. //    It suffers from 2 problems: a) Independence
  46. //    isnt great, that is the (n+1)th number is
  47. //    somewhat related to the preceding one, unlike
  48. //    flipping a coin where knowing the past outcomes
  49. //    dont help to predict the next result.  b)
  50. //    Taking parts of a LCG generated number can be
  51. //    quite non-random: for example, looking at only
  52. //    the least significant byte gives a permuted
  53. //    8-bit counter (that has a period length of only
  54. //    256).  The advantage of an LCA is that it is
  55. //    perfectly uniform when run for the entire period
  56. //    length (and very uniform for smaller sequences
  57. //    too, if the parameters are chosen carefully).
  58. //    
  59. //    ACGs have extremly long period lengths and
  60. //    provide good independence.  Unfortunately,
  61. //    uniformity isnt not too great. Furthermore, I
  62. //    didnt find any theoretically analysis of ACGs
  63. //    that addresses uniformity.
  64. //    
  65. //    The RNG given below will return numbers
  66. //    generated by an LCA that are permuted under
  67. //    control of a ACG. 2 permutations take place: the
  68. //    4 bytes of one LCG generated number are
  69. //    subjected to one of 16 permutations selected by
  70. //    4 bits of the ACG. The permutation a such that
  71. //    byte of the result may come from each byte of
  72. //    the LCG number. This effectively destroys the
  73. //    structure within a word. Finally, the sequence
  74. //    of such numbers is permuted within a range of
  75. //    256 numbers. This greatly improves independence.
  76. //    
  77. //
  78. //  Algorithm M as describes in Knuths "Art of Computer Programming",
  79. //    Vol 2. 1969
  80. //  is used with a linear congruential generator (to get a good uniform
  81. //  distribution) that is permuted with a Fibonacci additive congruential
  82. //  generator to get good independence.
  83. //
  84. //  Bit, byte, and word distributions were extensively tested and pass
  85. //  Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity
  86. //  assumption holds with probability > 0.999)
  87. //
  88. //  Run-up tests for on 7E8 numbers confirm independence with
  89. //  probability > 0.97.
  90. //
  91. //  Plotting random points in 2d reveals no apparent structure.
  92. //
  93. //  Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i),
  94. //    i=1..512)
  95. //  results in no obvious structure (A(i) ~ const).
  96. //
  97. //  Except for speed and memory requirements, this generator outperforms
  98. //  random() for all tests. (random() scored rather low on uniformity tests,
  99. //  while independence test differences were less dramatic).
  100. //
  101. //  AGN would like to..
  102. //  thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help.
  103. //
  104. //  And I would (DGC) would like to thank Donald Kunth for AGN for letting me
  105. //  use his extensions in this implementation.
  106. //
  107.  
  108. //
  109. //    Part of the table on page 28 of Knuth, vol II. This allows us
  110. //    to adjust the size of the table at the expense of shorter sequences.
  111. //
  112.  
  113. static randomStateTable[][3] = {
  114. {3,7,16}, {4,9, 32}, {3,10, 32}, {1,11, 32}, {1,15,64}, {3,17,128},
  115. {7,18,128}, {3,20,128}, {2,21, 128}, {1,22, 128}, {5,23, 128}, {3,25, 128},
  116. {2,29, 128}, {3,31, 128}, {13,33, 256}, {2,35, 256}, {11,36, 256},
  117. {14,39,256}, {3,41,256}, {9,49,256}, {3,52,256}, {24,55,256}, {7,57, 256},
  118. {19,58,256}, {38,89,512}, {17,95,512}, {6,97,512}, {11,98,512}, {-1,-1,-1} };
  119.  
  120. //
  121. // spatial permutation table
  122. //    RANDOM_PERM_SIZE must be a power of two
  123. //
  124.  
  125. #define RANDOM_PERM_SIZE 64
  126. unsigned long randomPermutations[RANDOM_PERM_SIZE] = {
  127. 0xffffffff, 0x00000000,  0x00000000,  0x00000000,  // 3210
  128. 0x0000ffff, 0x00ff0000,  0x00000000,  0xff000000,  // 2310
  129. 0xff0000ff, 0x0000ff00,  0x00000000,  0x00ff0000,  // 3120
  130. 0x00ff00ff, 0x00000000,  0xff00ff00,  0x00000000,  // 1230
  131.  
  132. 0xffff0000, 0x000000ff,  0x00000000,  0x0000ff00,  // 3201
  133. 0x00000000, 0x00ff00ff,  0x00000000,  0xff00ff00,  // 2301
  134. 0xff000000, 0x00000000,  0x000000ff,  0x00ffff00,  // 3102
  135. 0x00000000, 0x00000000,  0x00000000,  0xffffffff,  // 2103
  136.  
  137. 0xff00ff00, 0x00000000,  0x00ff00ff,  0x00000000,  // 3012
  138. 0x0000ff00, 0x00000000,  0x00ff0000,  0xff0000ff,  // 2013
  139. 0x00000000, 0x00000000,  0xffffffff,  0x00000000,  // 1032
  140. 0x00000000, 0x0000ff00,  0xffff0000,  0x000000ff,  // 1023
  141.  
  142. 0x00000000, 0xffffffff,  0x00000000,  0x00000000,  // 0321
  143. 0x00ffff00, 0xff000000,  0x00000000,  0x000000ff,  // 0213
  144. 0x00000000, 0xff000000,  0x0000ffff,  0x00ff0000,  // 0132
  145. 0x00000000, 0xff00ff00,  0x00000000,  0x00ff00ff   // 0123
  146. };
  147.  
  148. //
  149. //    SEED_TABLE_SIZE must be a power of 2
  150. //
  151. #define SEED_TABLE_SIZE 32
  152. static unsigned long seedTable[SEED_TABLE_SIZE] = {
  153. 0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b,
  154. 0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf,
  155. 0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706,
  156. 0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10,
  157. 0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a,
  158. 0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32,
  159. 0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f,
  160. 0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf
  161. };
  162.  
  163. //
  164. //    The LCG used to scramble the ACG
  165. //
  166. //
  167. // LC-parameter selection follows recommendations in 
  168. // "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi.
  169. //
  170. // LC_A = 251^2, ~= sqrt(2^32) = 66049
  171. // LC_C = result of a long trial & error series = 3907864577
  172. //
  173.  
  174. static const unsigned long LC_A = 66049;
  175. static const unsigned long LC_C = 3907864577;
  176. static inline unsigned long LCG(unsigned long x)
  177. {
  178.     return( x * LC_A + LC_C );
  179. }
  180.  
  181.  
  182. ACG::ACG(unsigned long seed, int size)
  183. {
  184.  
  185.     initialSeed = seed;
  186.     
  187.     //
  188.     //    Determine the size of the state table
  189.     //
  190.     
  191.     for (register int l = 0;
  192.      randomStateTable[l][0] != -1 && randomStateTable[l][1] < size;
  193.      l++);
  194.     
  195.     if (randomStateTable[l][1] == -1) {
  196.     l--;
  197.     }
  198.  
  199.     initialTableEntry = l;
  200.     
  201.     stateSize = randomStateTable[ initialTableEntry ][ 1 ];
  202.     auxSize = randomStateTable[ initialTableEntry ][ 2 ];
  203.     
  204.     //
  205.     //    Allocate the state table & the auxillary table in a single malloc
  206.     //
  207.     
  208.     state = new unsigned long[stateSize + auxSize];
  209.     auxState = &state[stateSize];
  210.  
  211.     reset();
  212. }
  213.  
  214. //
  215. //    Initialize the state
  216. //
  217. void
  218. ACG::reset()
  219. {
  220.     register unsigned long u;
  221.  
  222.     if (initialSeed < SEED_TABLE_SIZE) {
  223.     u = seedTable[ initialSeed ];
  224.     } else {
  225.     u = initialSeed ^ seedTable[ initialSeed & (SEED_TABLE_SIZE-1) ];
  226.     }
  227.  
  228.  
  229.     j = randomStateTable[ initialTableEntry ][ 0 ] - 1;
  230.     k = randomStateTable[ initialTableEntry ][ 1 ] - 1;
  231.  
  232.     register int i;
  233.     for(i = 0; i < stateSize; i++) {
  234.     state[i] = u = LCG(u);
  235.     }
  236.     
  237.     for (i = 0; i < auxSize; i++) {
  238.     auxState[i] = u = LCG(u);
  239.     }
  240.     
  241.     k = u % stateSize;
  242.     int tailBehind = (stateSize - randomStateTable[ initialTableEntry ][ 0 ]);
  243.     j = k - tailBehind;
  244.     if (j < 0) {
  245.     j += stateSize;
  246.     }
  247.     
  248.     lcgRecurr = u;
  249.     
  250.     assert(sizeof(double) == 2 * sizeof(long));
  251. }
  252.  
  253. ACG::~ACG()
  254. {
  255.     if (state) delete state;
  256.     state = 0;
  257.     // don't delete auxState, it's really an alias for state.
  258. }
  259.  
  260. //
  261. //    Returns 32 bits of random information.
  262. //
  263.  
  264. unsigned long ACG::asLong()
  265. {
  266.     unsigned long result = state[k] + state[j];
  267.     state[k] = result;
  268.     j = (j <= 0) ? (stateSize-1) : (j-1);
  269.     k = (k <= 0) ? (stateSize-1) : (k-1);
  270.     
  271.     short int auxIndex = (result >> 24) & (auxSize - 1);
  272.     register unsigned long auxACG = auxState[auxIndex];
  273.     auxState[auxIndex] = lcgRecurr = LCG(lcgRecurr);
  274.     
  275.     //
  276.     // 3c is a magic number. We are doing four masks here, so we
  277.     // do not want to run off the end of the permutation table.
  278.     // This insures that we have always got four entries left.
  279.     //
  280.     register unsigned long *perm = & randomPermutations[result & 0x3c];
  281.     
  282.     result =  *(perm++) & auxACG;
  283.     result |= *(perm++) & ((auxACG << 24)
  284.                | ((auxACG >> 8)& 0xffffff));
  285.     result |= *(perm++) & ((auxACG << 16)
  286.                | ((auxACG >> 16) & 0xffff));
  287.     result |= *(perm++) & ((auxACG <<  8)
  288.                | ((auxACG >> 24) &   0xff));
  289.     
  290.     return(result);
  291. }
  292.